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The limits to the error due to truncation of 
the numeric integration of the one-sided 
Laplace transform of a relaxation function 
in the time domain into its equivalent fre- 
quency domain are established. Separate 
results are given for large and small o). 
These results show that, for a given o), only 
a restricted range of time samples is needed 
to perform the computation to a given accu- 
racy. These results are then combined with 
a known error estimate for integration by 
cubic splines to give a good estimate for 
the number of points needed to perform the 



computation to a given accuracy. For a 
given data window between fi and f2, the 
computation time is shown to be propor- 
tional to ln(fi/f2). 
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1. Introduction 

The transformation of relaxation data recorded in the 
time-domain into their equivalent in the frequency- 
domain has been a difficult problem. The data can cover 
a wide range of times with a correspondingly wide 
range of frequencies. Both frequency and time are typ- 
ically measured on logarithmic scales that can cover 
many decades. What would be desirable is a method for 
computing the transform to a level of accuracy consis- 
tent with the original time data. Also, the time and space 
requirements should be such that transform can be 
performed as the data are being acquired so that the 
frequency domain data can be followed during the 
course of a long measurement in the time domain. 

In a previous paper [1], it was demonstrated that the 
use of a cubic spline as an integrating tool provided a 
rapidly convergent method as a function of sampling 
density for obtaining the Laplace transform of a relax- 
ation function measured as a function of time. The data 
were assumed to be sampled uniformly on a logarithmic 
time scale and that the results were displayed on the 
equivalent logarithmic frequency scale. The method 
was numerically stable and quite efficient, being order 



0[ln{t2/tif] in computation, where ti and ti are the 
beginning and ending times of the measurement. 

Some problems remained in the method. Unlike a 
conventional FFT type transform, the logarithmic spac- 
ing requires the integration of sine and cosine functions 
that are not harmonically related. Also, the integrals 
over powers times the sines and cosines must be evalu- 
ated very accurately for small arguments where the 
analytic expressions break down numerically. These 
conditions require either an extensive computation per 
point, slowing down the transform, or else the use of a 
look-up table that can get very large. An examination of 
the quantities involved in the computation suggested 
that too much computation was being done. 

This paper demonstrates that it is possible to restrict 
the computation to a predetermined diagonal band 
transformation matrix such that the computation be- 
comes 0[ln{t2/ti)] and thus far more efficient. As a 
consequence, the transformation need be defined for at 
most a normalized frequency decade. This drastically 
reduces the size of the look-up table so that it can be 
readily precomputed prior to the main computation, or 
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stored as a table without requiring too much space. 
Furthermore, the computation becomes almost indepen- 
dent of the size of the data window. 

The banding criteria will be developed in terms of a 
desired tolerance. This will be done separately for both 
small and large (at. When combined with an estimate for 
the point density needed to sample a relaxation function 
with a spline, these results lead to an explicit estimate 
for the number of points in the band. 

It should be noted that these results do not address the 
problem of continuing the integral outside the data lim- 
its. This problem is independent of the computation con- 
sidered here and has been addressed previously [1], [2]. 
However, the results presented here do address the ques- 
tion as to how much the missing data outside the mea- 
suring window affects the result for any computed fre- 
quency. 
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For a relaxation function that can be represented as a 
distribution of exponentials, all derivatives exist. Fur- 
thermore, if the samples used to determine S are uni- 
formly spanned in In t, then the mesh spacing h uni- 
formly approaches zero as the sampling density is 
increased. This leads to the result that 
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where cr is a tolerance and h = \n (tj/tj.i) [4]. For an 
exponential, the constant of proportionality for base 10 
logorithms can be taken as unity direct computation [1], 
so that 



2. Convergence Criteria 

Cubic splines have been intensively studied for their 
convergence properties. These results can be applied 
directly to the transform of a relaxation function. For the 
problem considered here, consider a relaxation function 
C{t) normalized to vary between and 1, measured in 
a data window /i < / < /j- One wishes to compute the 
Laplace transform of C{t+) numerically: 



(T^[lg{tj/tj.,)T 
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In general, a relaxation function is not a simple expo- 
nential. It can, however, be represented as a distribution 
of exponentials [5]. For any given exponential, the error 
term represented in Eq. (6) is a strongly peaked function 
of (w in 1 < u)T< 10. Therefore, for a general relaxation 
function and for a given w, only that part of the distribu- 
tion with values of wt near the peak in the error term 
can contribute as much as Eq. (6). Therefore, the error 
for a relaxation function characterized by a broad distri- 
bution of relaxation times must be less than that given 
by Eq. (6). 



Real(5) = 0+. 



(1) 



3. Sliort Time Limit 



In the data window this becomes 



C* = 



e"'"' C At, 



(2) 



where C is the first derivative of C with respect to t. 

Let Sih) be the spline passing through Cit) using a 
mesh spacing measured by h . Then the integral can be 
written as 



: {e"' S'{h)At+ I 



C*(w) = e"" S'{h) At + e"" (C - S'{h)At, (3) 



where the second term represents the error from using 
the spline to represent the relaxation function. Since 
|e"""'| = 1, the second term can be bounded by 



Since a cubic spline is a piecewise polynomial, once 
the spline is fitted to the data, one must evaluate the 
integrals of the form 



x" sin(x) dx, 

x" sin(x) dx, 
n = 0, 1,2. 



(7) 



Only powers to the second degree are needed as the 
integration is carried out on the derivative of the fitted 
spline. As noted before [1], the best way to evaluate 
these integrals for small x = wt is to use the Taylor 
expansion of the integrals. If only the leading term of 
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the Taylor expansions are kept, except for the intgral over 
cos(x), the relative error R is given by 



« < y, x< 1. 



(8) 



For a full scale of unity 



■re" 

T Jo 



dt= 1 - e" 



This shows that the explicit integration can be cut off at 

X = la-'"- . (9) 

To this integration, a running sum of the derivative C 
must be added to the cosine integral which represents 
the real component of transform. 

From Eq. (6), the number of points required per 
decade is 



1 A""* 



(10) 



This then gives the number of data points A'l for x < 1 

as 



^^ = 2^H2b, 



(11) 



For the data corresponding to x smaller than the cut- 
off, one has simply 

. mo I \ 

C*(a;) = Q,«,) - C, - y 2 Q' (?,' - ^-i j , (12) 



where }(NC) is the highest index corresponding to 
wt < 2o""^. These terms can easily be kept as a running 
sum. 

4. Long Time Limit 

Whereas the short time data can be kept as a running 
sum, the long time data beyond a cutoff can be simply 
disregarded. To show this, consider the long time part of 
the transform 



•'ti 
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This integral is easily evaluated to give 
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This gives 
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The maximum value is give in the limit of y small, 
T » ti to give 



L., 
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Thus, for a tolerance of a, cot2 = llcr, one obtains for 
a cutoff from Eq. (18). 
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o- h'- 
This gives the number of data points A'h for x > 1 as 
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(20) 



5. Final Results 

The results given above can be simply combined to 
give the total width over which the data need to be 
integrated. This is give from Eqs. (1 1) and (20) as 



N-- 



T 
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This expression can be readily evaluated for a = 10'^, 
10"^, and 10'^ as 25, 59 and 131 data points respectively. 
For a given precision, the frequency spacing should 
correspond to the sampling density, which is a constant 
for a given precision. Therefore, the number of frequen- 
cies calculated is proportional to ln{t2/ti). Also, as the 
computation is a fixed number of operations for a given 
frequency, the total computation must be 0[ln{t2/ti)]. 
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This makes the algorithm only slowly dependent on the 
range of the data and, for wide frequency ranges, very 
efficient. 

In a computation on actual data, since the integrals 
can be defined in terms of cot, if both w and t are equally 
spaced in log w and log t, then only a single row of 
values need to be stored, since it can be shifted along the 
data. However, this typically results in awkward fre- 
quency values and convenient, rounded values are more 
desirable. In this case, one only has to use a table large 
enough to correspond to a normalized decade of fre- 
quencies. This is still a great savings in time and storage. 

Where the transformation band crosses the edges of 
the data at ti and t2 mark the points at which behavior 
outside the measurement window in the time domain 
affect the desired results in the frequency domain. For 
the parts of the desired frequency domain where this 
occurs, one must add the continuation mentioned previ- 
ously to minimize a possibly large error [2]. 

In the time-domain instrument in which this 
algorithm was embedded [2], experimental limitations 
prevented a true logarithmic spacing for the first few 
time samples as the sampling rate was not dense 
enough. This limitation was easily overcome by a sub- 
sidiary calculation over the first few data points, until 
the logarithmic spacing could be properly used, and then 
picking up the integral table at the appropriate point. 
While this did introduce some complexity to the pro- 
gram, it had little effect on the running time of the 
algorithm. For a spline taken at ten points per decade 
and for the transformation taken for cot from 10"^ to 10^, 
the computation was only a few times longer than the 
spline fitting, which used an over relaxation method that 
also was nearly 0[ln{t2/t,)] in computation. This band 
was a little broader than the 3X10"^ to lO'* band needed 
to match the maximum error of 1 0'' to minimize round- 
off and truncation error for the instrument. 

The algorithm has also been compiled into a 
FORTRAN subroutine' to carry out the numerical 
transform for a supplied real function. In this case, the 
time series data window is set by the required frequency 
window and tolerance. Also, the integration window has 
been made symmetric about the normalized frequency 
for clarity in the code so that running sums are not 
required and only the value of the transformed function 
at the beginning of the summation is needed for initial- 
ization. In this routine, the tolerance is a passed parame- 
ter and the routine determines the data window required 
to perform the transform. 

For illustration, the time function 



with known transform. 



C(t) 



■y(l3,t) 



(22) 



C*(w) 



1 



(l+ioJTf 



(23) 



the Cole-Davidson Function [7], where y{a,x) and F 
are the incomplete gamma and gamma functions respec- 
tively [8], was numerically transformed and compared 
with the analytic transform. For a requested tolerance of 
10'^, the maximum computed error was half that. The 
results are shown in Fig. 1. 
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Fig. 1. The computational error for a specified tolerance of 10 ^ for 
the Cole-Davidson relaxation function with an exponent of 0.5. 
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' The source code is available upon request from the author. 
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